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We present a method for the direct calculation of the spin stiffness by means of the coupled 
cluster method. For the spin-half Heisenberg antiferromagnet on the square, the triangular and the 
cubic lattices we calculate the stiffness in high orders of approximation. For the square and the 
cubic lattices our results are in very good agreement with the best results available in the literature. 
For the triangular lattice our result is more precise than any other result obtained so far by other 
approximate method. 
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I. INTRODUCTION 

The study of quantum magnetism has attracted much 
experimental and theoretical attention over many years, 
for an overview, see Ref. The spin stiffness p s con- 
stitutes, together with the spin-wave velocity, a funda- 
mental parameter that determines the low-energy dy- 
namics of magnetic systems^*^ In particular, in two- 
dimensional quantum antiferromagnets, where magneti- 
cally ordered as well as quantum disordered ground-state 
phases are observed, the ground-state stiffness measures 
the distance of the ground state from criticality^ and can 
be used, in addition to the sublattice magnetization M, 
to test the existence or absence of magnetic long-range 
order (LRO). 

Over the last 15 years in a series of pa- 
pers several methods like series expansion*^ spin- 
wave theory^i2i2iiSiii quantum Monte Carlo^ exact 
diagonalizationj&i&ii Schwinger-boson approachfi£ii£ii! 
and renormalization group theory 1 - have been used to 
calculate the spin stiffness of the spin-half Heisenberg 
antiferromagnet (HAFM) on the square, the triangular 
and the cubic lattices. However, results for the trian- 
gular lattice seem to be less precise than those for the 
square lattice due to strong frustration. Published val- 
ues therefore show significant variability. 

The spin stiffness p s measures the increase in the 
amount of energy when we rotate the order parameter of 
a magnetically long-range ordered system along a given 
direction by a small angle per unit length, i.e. 

N N 2 w 

where E{&) is the ground-state energy as a function of 
the imposed twist, and N is the number of sites. In the 
thermodynamic limit, a positive value of p s means that 
there is LRO in the system, while a value of zero reveals 
that there is no LRO(i£ 

In this paper we present a new method to calcu- 
late the spin stiffness for the quantum-spin HAFM us- 



ing the coupled cluster approachi22i2ii2£ The coupled 
cluster approach is a powerful and universal tool in 
quantum many-body physics which has been applied in 
various fields like nuclear physics, quantum chemistry, 
strongly correlated electrons etc.i^2i2L2^ More recently 
the coupled cluster method (CCM) has been applied 
to q uantum spin sy stems with much success, see, e.g., 
Refs. I23l24l25l26l2ll2l2l3cll3l3l3l3l In the field of 
magnetism an important advantage of this approach is 
its applicability to strongly frustrated quantum spin sys- 
tems in any dimension, where some other methods, such 
as, e.g., the quantum Monte Carlo method fail. There- 
fore the method to calculate the spin stiffness described 
in this paper is quite generally applicable to spin systems 
also with non-collinear ground states. 

To demonstrate the potential of the presented method 
we calculate the spin stiffness for the spin-i HAFM with 
nearest-neighbor interaction on the cubic, the square, 
and on the triangular lattices and compare our results 
with available data in the literature. While for the square 
and the cubic lattices accurate high order spin-wave re- 
sults are available which can be used to estimate the ac- 
curacy of the CCM results, the known results for the 
frustrated HAFM on the triangular lattice with a non- 
collinear ground state seem to be less reliable, since the 
used methods are less accurate. We argue that our result 
for the stiffness of the HAFM on the triangular lattice ob- 
tained by CCM in high order of approximation is better 
than the so far available results. We mention that some 
preliminary results for the spin stiffness of the so-called 
J — J' model using the CCM can be found in Ref. |32- 



II. THE METHOD 



The model we consider is the spin-half HAFM 



(2) 
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In J3J the sum runs over all pairs of nearest neighbors 
We now set J = 1 henceforth. 

We start with a brief illustration of the main features 
of the CCM. For a general overview on the CCM the 
interested reader is referred, e.g., to Ref. I22I and for de- 
tails of the CCM computational algorithm for quantum 
spin systems (with spin quantum number s = 1/2) to 
Refs. l26l28l29l The starting point for a CCM calculation 
is the choice of a normalized model or reference state |$), 
together with a set of mutually commuting multispin cre- 
ation operators C+ which are defined over a complete set 
of many-body configurations I. The operators C'i are the 
multispin destruction operators and are defined to be the 
Hermitian adjoints of the C+ . We choose {|$};C/} in 
such a way that we have ($1(7/ = = C/|$), V/ ^ 0. 
Note that the CCM formalism corresponds to the ther- 
modynamic limit N — > 00. 

For spin systems, an appropriate choice for the CCM 
model state |$) is often a classical spin state, in which 
the most general situation is one in which each spin can 
point in an arbitrary direction. 

We then perform a local coordinate transformation 
such that all spins are aligned in negative z-direction in 
the new coordinate frameiS&M As a result we have 

I*) = I • •• Hi ■■ •>; Cj = s+, a+8+, S+S+S+, 

(3) 

(where the indices i,j,k,... denote arbitrary lattice 
sites) for the model state and the multispin creation op- 
erators which now consist of spin-raising operators only. 
In the new coordinate system the Hamiltonian reads 2 ^ 




with tp being the angle between the two spins, and s = 
s x ± is v are the spin-raising and spin-lowering operators. 
According to Fig.^we have e.g. for the twisted Neel state 
on the square lattice tp = tt for nearest-neighbors along 
the y direction but tp = n + 9 along the x direction and 
for the twisted 120° Neel state on the triangular lattice 
we have tp — 2ir/3 + 9/2 for nearest neighbors along the 
^e x + e y direction but tp = Att/3 + 9 along the x direction. 

The CCM parameterizations of the ket and bra ground 
states are given by 

H\9) = E\$) ; (V\H = E(&\ ; 

|*} = e s |$) ; S = £><7+ i 

(*| = <$|£e- s ; S=1 + ^2S X CY. (5) 

IjiO 

The correlation operators S and S contain the correlation 
coefficients Si and Si that we must determine. Using 



the Schrodinger equation, H\^} = E\ty), we can now 
write the ground-state energy as E — ($>\e~ s He s \<S>) 
and the sublattice magnetization is given by M = 
— 1/NJ2^ (^\sf where sf is expressed in the trans- 
formed coordinate system. To find the ket-state and bra- 
state correlation coefficients we require that the expecta- 
tion value H = (iS?\H\ty) is a minimum with respect to 
the bra-state and ket-state correlation coefficients, such 
that the CCM ket- and bra-state equations are given by 

(<5>\Cjre- s He s \<S>) = 11^0 (6) 
{<l>\Se- s [H, C+]e s |$) = VJ + 0. (7) 

The problem of determining the CCM equations now be- 
comes a pattern-matching exercise of the {Cj} to the 
terms in e~ He s in Eq. JHJ). 

The CCM formalism is exact if we take into account all 
possible multispin configurations in the correlation oper- 
ators S and S, This is, however, generally not possible for 
most quantum many-body models including those stud- 
ied here. We must therefore use the most common ap- 
proximation scheme to truncate the expansion of S and S 
in the Eqs. © and J7J), namely the LSUBn scheme, where 
we include only n or fewer correlated spins in all config- 
urations (or lattice animals in the language of graph the- 
ory) which span a range of no more than n adjace nt ( con- 
tiguous) lattice sites (for more details see Refs. I2l2ll2l . 

The spin stiffness considered in this paper is the stiff- 
ness of the Neel order parameter (sublattice magnetiza- 
tion). Hence the corresponding model state |<E>) is the 
Neel state. This is the ordinary collinear two-sublattice 
Neel state for the square and the cubic lattices. The 
model state is a noncollinear 120° three-sublattice Neel 
state for the triangular lattice. Note that for the collinear 
Neel state only LSUBn approximations with even n are 
relevant ^SiSS In order to calculate the spin stiffness di- 
rectly using Eq. (JIJ we must modify the model (Neel) 
state by introducing an appropriate twist 9, see Fig. ^ 
Thus the ket-state correlation coefficients Si (after solv- 
ing the CCM equations ©) depend on 9 and hence the 
ground-state energy E is also dependent on 9. Note that 
our numerical code for the CCM-LSUBrt approximation 
allows us to calculate E(9) with very high precision of 
about 14 digits. First we have checked numerically that 
the ground-state energy calculated in LSUBn approxima- 
tion does indeed fulfill the relation Q with high precision 
for 9 < 0.01. The stiffness now can easily be calculated 
using numerical differentiation of E(9) which was done 
using a three-point formula with 9 = — 1 CP 4 , 0, +10 -4 . 

Since the LSUBn approximation becomes exact for 
n — > 00, it is useful to extrapolate the 'raw' LSUBn re- 
sults to the limit n — > 00. Although we do not know the 
exact scaling of the LSUBn results, there is some empir- 
ical experienc o 26 ' 28 i 29 how the ground-state energy and 
the order parameter for antiferromagnetic spin models 
scale with n. Based on this experience we have tested sev- 
eral fitting functions for the stiffness and we have found 
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TABLE I: Spin stiffness p s for the spin-half Heisenberg anti- 
ferromagnet on the square lattice calculated by various CCM- 
LSUBn approximations and the result of the n — » co extrap- 
olation using LSUBn with n = 4, 6, 8. 



LSUBn 


number eqs. 


Ps 


2 


3 


0.2574 


4 


40 


0.2310 


6 


828 


0.2176 


8 


21124 


0.2097 


extrapol. 




0.1812 



FIG. 1: Illustration of the twisted Neel state (a: square lat- 
tice; b: triangular lattice). The twist is introduced along rows 
in x direction. The angles at the lattice sites indicated the 
twist of the spins with respect to the corresponding Neel state. 



the best extrapolation is obtained by the fitting function 

1 1 



a = ao + a%- + a 2 — x. 

n n z 



(8) 



This law is k n0 wn^2i22i22i2ii22i2i to provide good extrap- 
olated results for the order parameter. We show this 
extrapolation in Fig. [21 



TABLE II: Collection of data for the spin stiffness p 3 for the 
spin-half Heisenberg antiferromagnet on the square lattice cal- 
culated by different methods. 



first-order spin-wave theory—*- 


0.191 


second-order spin-wave theory— i- 


0.181 


third-order spin-wave theory- 


0.175 


series expansion^ 


0.182 


exact diagonalization— 


0.183 


quantum Monte Carlo— 


0.199 


Schwinger- boson approach 1— 


0.176 


Schwinger-boson approach II— 


0.153 


CCM 


0.181 
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FIG. 2: Extrapolation of the CCM-LSUBn results for the 
stiffness. The points represent the CCM-LUBn results and 
the lines correspond to the function (|HJ fitted to this data 
points. 



III. RESULTS 

Let us start with the results for the square lattice. Ex- 
ploiting the lattice symmetries we are able to perform 
calculations up to LSUB8, where for the twisted state 
21124 ket equations © have to be solved. The results 



for the stiffness are given in Tab. [A Using LSUBn with 
n = 2,4,6,8 the extrapolated result is p s = 0.1831. As 
known from the sublattice magnetization even better re- 
sults can be obtained by excluding the LSUB2 data. In- 
deed the extrapolation using the LSUB4, LSUB6, LSUB8 
data yields p s — 0.1812. Note that the correspond- 
ing extrapolated value for the sublattice magnetization 29 
M = 0.3114 is in good agreement with other results^ 9 *^. 
A certain estimate of the accuracy can be obtained by an 
extrapolation using LSUB2, LSUB4, LSUB6, only, which 
yields p s = 0.1839. We compare our results for p s with 
some data obtained by other methods in Tab. |H| Ob- 
viously, there is a significant variance in the data. In 
particular, the value obtained by quantum Monte Carlo 
seems to be surprisingly large. However, this might be 
connected with the fact, that in Ref. [lj the stiffness 
was not determined directly, but via the temperature de- 
pendence of the correlation length which may lead to 
larger uncertainty. We think, that the high-order spin- 
wave theory^ is the most systematic approach, since one 
can see how the stiffness changes with increasing order 
of approximation. Assuming the third-order order spin- 
wave results as a benchmark we find that our CCM result 
deviates by about 3%. 

For the triangular lattice the twist we consider (see 
Fig. ^) corresponds to the in-plane spin stiffness. Due 
to the noncollinear structure of the three-sublattice Neel 
state also LSUBn approximations with odd n appear. 
Furthermore the number of ket equations in a certain 
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TABLE III: In-plane spin stiffness p s for the spin-half Heisen- 
berg antiferromagnet on the triangular lattice calculated by 
various CCM-LSUBn approximations and the result of the 
n — > oo extrapolation using LSUBn with n = 2,3, 4, 5, 6, 7. 



LSUBn 


number eqs. 


Ps 


2 


3 


0.1188 


3 


14 


0.1075 


4 


67 


0.0975 


5 


370 


0.0924 


6 


2133 


0.0869 


7 


12878 


0.0824 


extrapol. 




0.0564 



TABLE IV: Collection of data for the spin stiffness p s for the 
spin-half Heisenberg antiferromagnet on the triangular lattice 
calculated by different methods. 



exact diagonalization— 


0.05 


first-order spin-wave theory-*- 


0.080 


Schwinger- boson approach— 


0.088 


CCM 


0.056 



level of approximation becomes larger then for the square 
lattice and as a results the highest level of approximation 
we are able to consider is LSUB7. The results for differ- 
ent LSUBn approximations are given in Tab. lIIll The ex- 
trapolation of the LSUBn data according to Eq. JHJ with 
n = 2, 4, 6 leads to p s = 0.0604 and with n = 2, 3, 4, 5, 6, 7 
to p s = 0.0564. Again the difference in the two values 
can be considered as a certain estimate of the accuracy. 
As a byproduct of our high-order calculation we can give 
here improved values for the sublattice magnetization M. 
So far results for M up to LSUB62£*2i are published. We 
can add M = 0.3152 (LSUB7) and M = 0.3018 (LSUB8). 
The corresponding extrapolated value using Eq. . (JSJ) and 
LSUBn with n = 2,3,4,5,6,7,8 is M = 0.2134, which 
is close to spin-wav o 7 i 36 and Green's function Monte 
Carlo^ 7 . results. The small values of the stiffness and 
the order parameter in comparison with the square lat- 
tice are attributed to the frustration leading to a non- 
collinear ground state and in combination with quantum 
fluctuations to a drastic weakening of magnetic order in 
the spin-half HAFM. 

We compare our results for p s with available results 
from literature, see Tab. IIVI Comparing the methods 
used to calculate p s for the square lattice (Tab. HT|) and for 
the triangular lattice (Tab. llVl) we see that the results for 



the triangular lattice are much less relia ble, since here the 
accuracy of the methods used in Refs. I7l8lll7l is limited. 
Assuming the same tendency as for the square lattice we 
can expect that the first-order spin-wave value for ps^& 
becomes smaller (and therefore closer to our CCM result) 
going to second- and third-order spin-wave theories. We 
believe that our result is indeed of higher accuracy than 

TABLE V: Spin stiffness p 3 for the spin-half Heisenberg anti- 
ferromagnet on the cubic lattice calculated by various CCM- 
LSUBn approximations and the result of the n — + oo extrap- 
olation using LSUBn with n — 2, 4, 6. 



LSUBn 


number eqs. 


Ps 


2 


4 


0.2527 


4 


106 


0.2416 


6 


5706 


0.2380 


extrapol. 




0.2312 



data for p s so far available. 

We present now our results for p s for the simple cubic 
lattice, see Tab. Here the highest level of approxima- 
tion we can consider is LSUB6. From Fig. [3 it becomes 
obvious, that there is only a weak dependence on the level 
of CCM approximation n. Therefore we expect that the 
extrapolation according to Eq. JSJ yielding p s = 0.2312 
is particular accurate. Indeed we find that our result is in 
very good agreement with the result obtained by second- 
order spin- wave theory 6 p s = 0.2343. Note that the l/s 
spin-wave expansion seems to converge very rapidly 6 - and 
therefore the second-order spin-wave theory is expected 
to yield a very precise result for p s . For the sublat- 
tice magnetization a corresponding extrapolation leads 
to M — 0.418122 coinciding to 1% with the high preci- 
sion third-order spin- wave result 



IV. SUMMARY 

In summary, we have presented a method for the direct 
calculation of the spin stiffness within the framework of 
the coupled cluster method. We obtain accurate values 
for the stiffness by applying this algorithm to high or- 
ders of LSUBn approximations for the spin-half isotropic 
Heisenberg antiferromagnet on various lattices with and 
without frustration. 
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